Nonlinear dielectric effect of dipolar fluids 
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The nonlinear dielectric effect for dipolar fluids is studied within the framework of the mean 
spherical approximation (MSA) of hard core dipolar Yukawa fluids. Based on earlier results for the 
electric field dependence of the polarization our analytical results show so-called normal saturation 
effects which are in good agreement with corresponding NVT ensemble Monte Carlo simulation data. 
The linear and the nonlinear dielectric permittivities obtained from MC simulations are determined 
from the fluctuations of the total dipole moment of the system in the absence of an applied electric 
field. We compare the MSA based theoretical results with the corresponding Langevin and Debye- 
Weiss behaviors. 



I. INTRODUCTION 



The studies of dielectric polarization in fluids [l| are 
based on the following relation between the polarization 
P and the electric field strength E inside the dielectric: 



AttP = (e E - 1)E, 



(1) 



where e E is the field dependent dielectric permittivity. 
The internal field is often called Maxwell field [H,0] which 
differs from the external field Eo applied to the dielectric 
medium. In the low-field limit the linear dielectric per- 
mittivity eo of an isotropic fluid is given by the ratio of 
the polarization to the internal field strength: 
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(2) 



E=0 



In strong electric fields the polarization response acquires 
in addition nonlinear contributions, so that the elec- 
tric permittivity turns into a nonlinear function of the 
Maxwell field: 



e E = eo + e 2 £ 2 + e 4 £ 4 



(3) 



For simple molecular liquids the coefficients of the power 
series decrease rapidly |3j . Thus in general dielectric ex- 
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pcrimcnts can be well described in terms of the linear 
dielectric permittivity. For liquids consisting of small 
molecules this scries can be limited to the second term. 
The corresponding nonlinear dielectric effect (NDE) is 
determined by the coefficient e 2 of the contribution ~ E 2 
in Eq. ©: 



62 = 47rA 
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(4) 



where A is the nonlinear dielectric permittivity which, 
upon substituting Eq. ([3]) into Eq. ([T]), is given by 



A 
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(5) 
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Although the nonlinear dielectric behavior of liquids 
has enjoyed a long lasting scientific interest 0, 0] accu- 
rate measurements of NDEs have become possible only 
recently due to the development of new techniques [3j], 
which in an applied electric field are capable to separate 
from the NDE, e.g., the Joule effect due to heating of 
the sample induced by its finite conductivity. Experi- 
mentally, A is determined from the small change of the 
dielectric permittivity (Ace) induced by static 0, Q or 
pulsed 0,1^, 0, S] strong electric fields and is detected 
by a weak radio-frequency probing field. In nonlinear 
optics the NDE can also provide useful information for 
laser induced molecular reorientations in isotropic and 
liquid crystalline phases 0]. The first NDE measure- 
ments were carried out by Herweg in diethyl ether Q 
yielding e 2 < 0. If e 2 has a sign opposite to that of 
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eo > one often speaks of normal saturation because in 
this case the first correction term to the linear behavior 
P^E for E^O is in line with the levelling off at large 
E. This effect is mainly connected to the ordering of the 
orientation of dipoles in strong electric fields. The nega- 
tive value of €2 corresponds to the negative value of the 
third order term of the power expansion of the Langevin 
function for small E. Strongly dipolar liquids, such as 
nitrobenzene, show anomalous (positive) dielectric satu- 
ration because in such systems the external electric field 
influences the formation of antiparallel pairs of dipolar 
nitrobenzene molecules 0. Nonlinear dielectric effects 
have become also very useful for analyzing intermolec- 
ular association 10j], conformational equilibria II ill . and 
critical phenomena in liquids and liquid mixtures [l2l. [7^ . 
Recently the NDE has also been used to study isotropic 
- mesophase transitions of various liquid crystals [14 1 . 



The theoretical models of nonlinear dielectric phenom- 
ena are based on classical electrostatics and statistical 
mechanics of liquids. There have been early attempts by 
Debye [IB] , Onsager [l6| , and Kirkwood to calculate 
the nonlinear dielectric permittivity on the basis of phc- 
nomenological theories of dielectric continua. Nice sum- 
maries of these nonlinear theories can be found in Rcfs. 
[l| and 18 1 . From a microscopic point of view Rasaiah et 
al. [3] and Martina and Stell [20] have proposed a statis- 
tical mechanics description for NDE and electrostriction 
on the basis of a quadratic hypernetted chain approxi- 
mation. Using molecular dynamics (MD) simulations the 
dielectric saturation of water has been studied by Alper 
and Levy [2l| , but they have not published any numerical 
value for the nonlinear dielectric permittivity. Yeh and 
Berkowitz 22| have obtained new MD simulation data for 
the electric field dependence of the dielectric permittivity 
of water and they found that e e decreases with increasing 
applied field strength in accordance with the expectation 
of a normal saturation effect. Their external field simu- 
lation data show good agreement with a phcnomcnolog- 
ical equation proposed by Booth . This equation has 
been used in the calculation of dielectric saturation of wa- 



ter in membrane protein channels [24(. Recently Fulton 
251 ] has compared the nonlinear dielectric permittivity of 
water obtained from simulation, theory, and experiment. 
He concluded that an upgraded approach by Booth [26 1 
renders the best agreement between the simulation and 
experimental data; but the extent of the agreement de- 
pends on the water model used for determining the cor- 
responding correlation functions. It is also shown that 



the calculated nonlinear dielectric permittivity strongly 
depends on the application of phcnomenological cavity 
and reaction field corrections. These findings underscore 
the ongoing interest in statistical mechanics analyses of 
nonlinear dielectric effects of dipolar liquids. 

Within the framework of density functional theory 
(DFT) and the mean spherical approximation (MSA) we 
have proposed an equation [271 ] for the magnetic field 
dependence of the magnetization of ferrofluids, which 
turned out to be successful in comparison with corre- 
sponding Monte Carlo (MC) simulation data. Translated 
into the synonymous electric language this means that we 
have obtained an analytical (implicit) equation for the 
electric field dependence of the polarization. Here, from 
our previous results [27| , we deduce a formula for the non- 
linear dielectric permittivity of dipolar fluids. Moreover, 
we compare our theoretical findings with canonical MC 
simulation data using the corresponding fluctuation for- 
mulae in the absence of external fields. The motivation 
for our study is to deduce an analytical theory of a re- 
alistic model for nonlinear phenomena in dipolar liquids. 
This theory has to be quantitatively correct in the sense 
that it withstands the comparison with MC simulation 
data and it has to be easily applicable for interpreting 
corresponding experimental data. We expect that our 
model calculations can be used to shed light on the role 
of nonlinear dielectric saturation for solvation effects of 
ions in dipolar solvents 



II. THEORY 



A. Microscopic model 



We study hard core dipolar Yukawa fluids which con- 
sist of spherical particles interacting via a hard core 
Yukawa (Y) potential characterized by parameters er, ey, 
and k: 



u Y {r\2) 



00 , 7-12 < a 

~e Y a{r 12 Y X exp[-K(ri 2 - a)\ , r 12 > a. 

(6) 

In addition there is a dipolar interaction due to point 
dipoles embedded at the centers of the particles: 



ur>(ri2, wi,w 2 ) = — o-£>(wi2, uji, uj 2 ), (7) 

?12 
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with the rotationally invariant function 



function of the DY fluid can be expressed as 



D(wi 2 , Wi, uj 2 ) = 3(fhi •fi 2 )(m 2 • r i2 ) - (mi • fh 2 ), (8) 

where particle 1 (2) is located at ri (r 2 ) and carries a 
dipolc moment of strength m with an orientation given 
by the unit vector mi(wi) (fh 2 (o; 2 )) with polar angles 
uji = (6»i, 0i) (w 2 = (6> 2 ,0 2 )); r 12 = r x - r 2 is the dif- 
ference vector between the centers of particle 1 and 2, 
12 = I r i2|, and ?i2 = T\i/ri2 is a unit vector with ori- 
entation lo\2 = (6*i 2 , 0i2). The hard core dipolar Yukawa 
interaction potential is defined by the sum of the afore- 
mentioned potentials as 

W£>y(ri2,wi,W2) = uy(ri 2 ) + Uu(r 12 , uji, w 2 ). (9) 



B. MSA for dipolar Yukawa fluids 

The MSA is defined by three equations relating the to- 
tal correlation function /i(ri 2) w 2 ) and the direct cor- 
relation function c(ri 2 , u>i, o; 2 ) as follows: 



4tt 



/i(ri2,wi,Lj 2 ) = c(ri 2 ,cji,cj 2 ) + 
dui 3 J d 3 r 3 h(r 13 , uj x , w 3 )c(r 23 , w 2 , w 3 ), 

/i(ri2,o;i,W2) = -1, r\2 < a, 



(10) 



(11) 



CDY (ri2, Wl, CJ 2 ) = cy(r i2 ) + 
CD(ri2)-D(Wl2,Wl,W2) + CA(ri 2 )A(cJi, w 2 ), (13) 



where 

A(wi,w 2 ) = mi-m 2 (14) 

is a rotationally invariant function. (We note that a sim- 
ilar equation is valid for the total correlation function 
h<DY (ri2) Wi, W2) of the DY fluid.) The radially symmet- 
ric function Cy(?"i 2 ) is the solution of the OZ equation 
with a simple hard-core Yukawa MSA closure and is given 
111 Refs. |35| while the functions CD(r 12 ) and CA(fi 2 ) 
are the solutions of two MSA integral equations (derived 
from the OZ equations) with the corresponding dipolar 
hard sphere MSA closure, depending therefore on the 
dipole moment of the molecules |3l| . In Eq. (fTB")) D and 
A are given by Eqs. ([8]) and (fT4]) respectively. The coef- 



ficients cd (7-12) and ca(^i 2 ) (together with /id ('"12) and 
/ia(^i 2 ) ) are independent of the Yukawa potential pa- 
rameters £y and k as well as of wi 2 , u>i, 0J2 and are given 
in Ref. [31|. The main feature of this solution is that 
it decomposes into contributions from the Yukawa po- 
tential and from the dipolar hard sphere potential with 
the latter ones factorizing into radial and angular depen- 
dences. In Ref. 33j it has been shown that, within the 
framework of MSA, the free energy Fdy of the DY fluids 
can be written as 



and 



c(ri 2 , U)\,u)2) = -/3u(ri2, r 12 >a. (12) 

Here and in the following j3 = 1/ksT, where &g is 
the Boltzmann constant and T is the temperature, and 
w(ri2, u>i, 0J2) is the pair potential characterizing the sys- 
tem. The number density p = N/V is given by the num- 
ber N of molecules in the system of volume V. Equation 
(fT0|) is the Ornstein-Zernike (OZ) relation, Eq. (|11|) is an 
exact relation for hard spheres, and Eq. (fT2|) is the clo- 
sure relation for the MSA (for details see Ref. [i3|). For 
dipolar hard sphere (DHS) fluids an analytical solution 
of the MSA was reported by Wertheim [31| . Later the 
MSA for hard core Yukawa fluids has also been solved 



analytically by Waisman 32|. Following the ideas of 
Wertheim, for hard core dipolar Yukawa fluids Henderson 
et al. 33[ 34| found an analytical solution in the frame- 
work of MSA. Within this theory the direct correlation 



Fdy — F, 



ID 



TpeX 771 

r DY - r ID 



*HS 



Fy X 



F d x hs, (15) 



where Fj d is the ideal gas free energy and Fff s , F Y X , and 



F D X HS are the excess free energies of hard sphere |3C 
hard core Yukawa [32| . and dipolar hard sphere fluids 
3l|, respectively. Within this approximation the dielec- 
tric constant of the DY fluid is given by_the formula due 
to Wertheim 



3l| and Henderson et al. [33j : 



where 



eo 



q{x) 



g(2g(2/)) 



(16) 



(17) 



is the reduced inverse compressibility function of the hard 
sphere fluid. The parameter £(y) stems from the DHS 
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MSA and is given by the implicit equation 

3y = q(20-q(-0, (18) 

where 

4-7T m 2 p 

v = l>]J (19) 

measures the reduced dipole strength. 

C. Field dependence of polarization 

The disadvantage of the MSA is that it can predict 
the polarization only for weak electric fields because in 
essence it is a linear response theory. In order to over- 
come this shortcoming we resort to density functional 
theory. Within this framework we first have to deter- 
mine the oricntational distribution function in an applied 
external field. For an inhomogencous and anisotropic 
dipolar fluid the one-particle distribution function p(r, uj) 
depends on the position r and the orientation lu of the 
particles. Accordingly the number density is given by 
p(r) = J duip(r,uj) and a(r, u) = p(r,uj)/p(r) is the 
oricntational distribution function. In a homogeneously 
polarized bulk phase (concerning influences from the 
sample shape see below) the number density p is spa- 
tially constant and the orientational distribution func- 
tion a(r, u) = a (8) depends only on the angle 6, which 
measures the orientation of the dipole of a single particle 
relative to the field direction. (Here we do not consider 
electrostriction, i.e., the dependence of p on E, because 
this contribution to Ae is one order of magnitude smaller 
than the corresponding contribution of the orientational 
ordering of dipoles in the presence of an applied electric 
field [36|). The orientational distribution function can 
be obtained by minimizing the following grand canonical 
variational functional: 



ft[p,{a(9)},T,p} 



= F ID [p,{a(6)},T} + 
F D * Y { Pl {a(9)},T}- 



trie field direction) so that the internal field strength E 
is equal to the strength of the externally applied electric 
field Eq (see Ref. [2[). The excess DY free energy func- 
tional for an anisotropic system is not known. However, 
it can be approximated by a functional Taylor series, ex- 
panded around a homogeneous isotropic reference system 
with bulk number density p. Neglecting all terms beyond 
second order, one has 

PF^ Y [p, {a(9)}, T] = pF%y{p, T) - 
y J (PrxdbJx j d 3 r 2 duj2^a(9 1 )Aa(02) X 

CDY(Tl2,Wl,U2), (21) 

where Aa(8) = a(0) — 1/(4tt) is the difference between 
the actual anisotropic and the isotropic orientational dis- 
tribution function and F D X Y (p? T) is the excess dipolar 
Yukawa free energy for an isotropic distribution given by 
the last three terms in Eq. (JT5J). Using the MSA direct 
correlation function of the DY fluid (Eq. I|13p) within this 
approximation the anisotropic excess free energy func- 
tional F D x Y [p, {a(9)},T] can be calculated analytically, 
which allows one to obtain a{9) by minimizing the grand 
canonical functional in Eq. (|2"D|) . From the orientational 
distribution function the polarization P (the direction of 
which coincides with the direction of the external electric 
field) can be obtained as 



P = p J duja{6)m cos 0. 



(22) 



Based on the orientational distribution function a{6) Eq. 
(|22p leads to the following polarization function: 



P = mpL (3mE + 3P 



(1 -«(-£(!/))) 



mp 



(23) 



p I d 3 rdu>a(9)(p + mE cost 



(20) 



where L(x) = cothx — 1/x is the well known Langevin 
function [l| . This is an implicit equation for the external 
field dependence of the polarization. The details of this 
calculation can be found in Ref. j2?| (where the equiv- 
alent magnetic language was adopted, i.e., the particles 
carry magnetic dipole moments and interact with an ap- 
plied magnetic field). 



where Fid and F D X Y are the ideal gas and the excess 
dipolar Yukawa free energy density functionals, respec- 
tively, p is the chemical potential, and E is the internal 
electric field in the sample. In this theory we assume that 
the volume of the fluid V has the shape of a macroscopic 
prolate rotational ellipsoid (elongated around the elec- 



D. Nonlinear dielectric effect 

In order to obtain the nonlinear dielectric permittivity 
A the third order derivative of the polarization with re- 
spect to the Maxwell field has to be calculated (see Eq. 



5 



([5])). To this end we introduce the following dimension- 
less form of Eq. (ggj) : 



where 



p(e) =L(e + Tp(e)) , 



p = P I {mp) 1 e = (3mE 



(24) 



(25) 



and 



d 3 p 



54 



e=0 



5 (3-r) 4 



I5g 4 (-£(y)) 



(32) 



Equations ([2|), (fl8|) . and (f30|) render back the expression 
in Eq. (fT6|) for the linear dielectric permittivity. For the 
nonlinear dielectric permittivity Eqs. and (f3"2")) yield 



arc the dimensionlcss polarization and electric field 
strength, respectively, and T = 3(1 — q{— £(y))) is a field 
independent parameter. From Eq. (f2"4")) one obtains 



dp dL(x) 
de dx 



(26) 



where a = e + Tp(e). The corresponding second order 
derivative is 



d 2 p _ d 2 L(x) 
de 2 dx 2 



i + r 



dp 
~de 



,dL(x) 



dx 



fp 
de^' 
(27) 



For the third order derivative one has 



d 3 p _ d 3 L(x) 
de 3 dx 3 



dp 

i + r-f 

de 



^dL{x) 



dx 



d 3 p 
de 3 



3F 



d 2 L(x) 



dx 2 



^ j,dp\ d 2 p 
de I de 2 



±5- (28) 



Since we consider thermodynamic states without sponta- 
neous polarization (see Ref. 27J) one has L(x = 0) = 
so that p(e = 0) = 0. With 



dL(x) 



dx 



x=0 



d 2 L(x) 



dx 2 
d 3 L(x) 



dx 3 



x=0 



one finds 



dp 
de 



e=0 



3-r U-ttv)Y 



fp 

de 2 



= o, 



(30) 



(31) 



e=0 



A 



™4 

m p 



1 



45(fc B T)3 <f{-t{y)) 

™2 



20n(k B T) 2 q^-^y))' 



(33) 



We note that for T = Eq. (|2"4"|) reduces to the Langevin 
equation for non-interacting dipoles, 



p(e) = L(e), 



(34) 



while for r = 3y one obtains the well known mean field 
Debye- Weiss polarization equation 



pie 



L(e + 3yp(e)). 



(35) 



Accordingly, for non-interacting dipoles (i.e., T = and 
(?(—£) = 1) one has for the linear and nonlinear dielectric 
permittivities 



A = 



^4 

m p 



45(fc B T)= 



20n(k B T) 



3y, 
y, 



(36) 



respectively, and within the framework of the Debye- 
Wciss theory (i.e., T = 3y and <?(—£) = 1 — y) 



eo = l + 3 



y 



i-y' 



4 

m p 



1 



45(fc B T) 3 (1 - y) 

™2 



20n(k B T) 2 (1 - y) 4 



(37) 



The latter equations show that not only the linear but 
also the nonlinear dielectric permittivity diverge for y = 
Airpm 2 1 '(9fcsT) — > 1, which for a fixed value m of the 
dipole moment provides the Debye- Weiss critical tem- 
perature k B T c = 4irpm 2 /9 of the isotropic liquid - fer- 
roelectric liquid second-order phase transition. In the 
following these theoretical predictions will be compared 
with corresponding Monte Carlo simulation results. 
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III. MONTE CARLO SIMULATIONS 



We have performed Monte Carlo simulations for DY 
fluids using the canonical NVT ensemble and applying 
Boltzmann sampling, periodic boundary conditions, and 
the minimum image convention [371 ]. In order to take 
into account the long-ranged character of the dipolar in- 
teraction the so-called reaction field (RF) method is used. 
According to this method we consider our system to be 
a macroscopic spherical sample composed of a number of 
replicas of the basic simulation cell embedded in a dielec- 
tric continuum with dielectric constant e RF . In this case 
in the spherical sample the internal (Maxwell) field E is 
not equal to the external applied field Eg, but (j^ 



E 



/ 2e RF 



1 



En 



(38) 



For such a system, using the third order Taylor series 
expansion of the polarization with respect to the external 
field, Kusalik [3911 has shown that 



Comparing Eq. (|42)) with Eq. ([4]) for the nonlinear di- 
electric permittivity we obtain 



A = 



90V 



(3<M 4 ) -5(M 2 ) 2 ) 



(44) 



In our NVT ensemble MC simulations e an d A are cal- 
culated from Eqs. (|4"3"| and (|4"4")h respectively. 

A spherical cutoff of the hard core DY pair potentials 
at half of the cubic box has been applied and long-ranged 
corrections (LRC) were taken into account 37j. In our 
simulations N — 256 particles have been used. We have 
not resorted to any finite-size scaling analysis for detect- 
ing the occurrence of the isotropic - anisotropic phase 
transitions. The simulations were started from a hep lat- 
tice configuration with randomly oriented dipolcs. After 
20.000 equilibration cycles, 2 x 10 6 - 4 x 10 6 production 
cycles were used. Statistical errors were calculated from 
the standard deviations of sub-averages containing 2x 10 5 
cycles. 



e B ~ e + 



2e RF + 1 



47T/3 3 

90V 



(3(M 4 ) -5(M 2 ) 2 ) 



where 



co = 



1 + 2e RF 



^^(M 2 )o 



3V 



l + 2e RF -^(M*) 



El (39) 



(40) 



In Eqs. (|3SJ) and |gDJ) M is the total dipole moment of 
the system of volume V, 



N 



M = ^m,, 



(41) 



and (M 2 )o, (M 4 )o are the ensemble averages of the cor- 
responding moments in zero external field. 

In our simulations we apply a conducting boundary 
condition which means e RF oo. In this limit Eqs. 
(|3"9")l . and ([4"0")) for the internal field dependent dielectric 
permittivity lead to 



4tt/? 3 
90V" 



(3(M 4 ) - 5(M 2 ) 2 ) E 1 



with 



eo = 1 



(42) 



(43) 



IV. RESULTS AND DISCUSSION 

In the following we shall use reduced quantities: T* = 
ksT/sY as the reduced temperature, p* = pa 3 as the re- 
duced density, m* = mj \Jey<y 3 as the reduced dipole mo- 
ment, and A* = Aey/cr 3 as the reduced nonlinear dielec- 
tric permittivity. Concerning the range of the Yukawa 
potential, all our results correspond to k = 1.8/ a. 

Figure 1 shows the linear and nonlinear dielectric per- 
mittivity as function of the reduced number density p*. 
For (m*) 2 = 0.5 and T* = 1 the predictions of our MSA 
theory (see Eqs. (fll))) - (fT5T) and (|3"3")) ) are compared with 
the Langcvin (Eq. i[36]l ) and Debye- Weiss (Eq. ([37]) ) 
approximations and our MC data. The Langevin theory 
predicts, in general, a linear dependence of eo and A on 
p. Figure 1 shows that the interparticle interaction en- 
hances eo and |A*| relative to the corresponding values 
of the Langcvin theory. The Dcbyc- Weiss theory over- 
estimates the MSA results for eo and |A*|. The critical 
density, at which within the Debye- Weiss theory eo and 
A diverge, is p c — 1.432. The behavior of eo and A shown 
in Fig. 1 for various approximations can be understood 
in terms of a low density expansion which amounts to an 
expansion in terms of y < 1 which is proportional to p 
(see Eq. (|T9| ). Within MSA the linear dielectric per- 
mittivity of the hard core dipolar Yukawa fluid has the 
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FIG. 1: Linear (a) and nonlinear (b) MSA dielectric proper- 
ties of hard core dipolar Yukawa fluids for (m*) 2 = 0.5 and 
T* = 1 in comparison with the corresponding Langevin and 
Debye-Weiss approximations. The symbols in both figures 
represent our Monte Carlo simulation data whereas the lines 
are obtained from the three corresponding theories. In (a) 
the size of the error bars is that of the symbols. Within the 
Debye-Weiss theory both e an d A* diverge at p* = 1.432. 
Note that according to Eqs. (I47|) and (|49[) A* vanishes lin- 
early for p* — > 0. 



expansion 34 1 



(e )MSA = 1 + 'iy + 3y 2 + -|</ 3 + 0(y 4 ). 

16 

From Eqs. flB]). ([15]). and (gSJ) one finds 

, t, = "V~ = 1 + v + JRV + 0{y 3 ), 
q(-£{y)) % 16 



(45) 



(46) 



A 



MSA 



20*(k B T)' 9»(-{(»)) 



The corresponding expansions within the Debye-Weiss 
approximation are 



(<*) 



DW 



1 + 3- 



y 



i-y 

1 + 3y + 3y 2 + 3y 3 + 0(y 4 ) = 

AK 

(e )MSA + YQ y 3 + 0(y 4 ) 



(48) 



and 



A 



DW 



y 



207r(fc B T) 2 

^MSA 



207r(fc B T) 2 (1 - y) 4 
(y + Ay 2 + I0y 3 ) + 0(y 4 ) 
3m 2 



167r(fc B T) 2 



y 3 + 0{y 4 ). 



(49) 



Equations (|48|) and (|49j) explain the trends of the various 
curves shown in Fig. 1. For the present choice of m* and 
T* the MC data and the MSA results for the linear dielec- 
tric permittivity agree very well. This agreement remains 
rather good also for the nonlinear dielectric permittivity. 
We note that the very good agreement between the MSA 
and the MC simulation data for the linear dielectric per- 
mittivity has also been reported in Ref. As shown 
in Fig. 2, the increase of the reduced temperature from 
T* = 1 to T* = 2 docs not change the character of the 
curves. The MSA results are in good agreement with the 
simulation data. With increasing temperature the dis- 
crepancy between the Debye-Weiss theory and the sim- 
ulation data is reduced because for higher temperatures 
the critical density, at which within the Dcbyc- Weiss the- 
ory eo and A diverge, is shifted to higher densities (here 
p* = 2.865; we note that this value of the density is phys- 
ically not accessible, because the corresponding packing 
fraction ?/ = irp* /6 would be larger than 1). The absolute 
values of eo and A arc considerably reduced upon rais- 
ing the temperature. Figure 3 shows that upon increas- 
ing the dipole moment from (m*) 2 = 0.5 to (m*) 2 = 1 
the agreement between the MSA and the Monte Carlo 
simulation data does not change significantly. However, 
the discrepancy between the Debye-Weiss theory and the 
simulation data has widened because due to the increase 
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of the dipole moment the critical density, at which within 
the Debye- Weiss theory eo and A diverge, is shifted to a 
lower value (/?* = 0.716). The Langcvin theory under- 
estimates both eo and |A*|. The increase of the dipole 
moment significantly increases eo and |A*|. 

Figure 4 shows that upon increasing the temperature 
from T* = 1 to T* = 2 at a fixed value (m*) 2 = 1 
of the dipole moment the agreement between the MSA 
and the simulation data is improved and the values of 
eo and |A*| are significantly reduced. The discrepancy 
between the Debye- Weiss theory and the simulation data 
is decreased because the critical density p* = 1.432, at 
which eo and A* diverge according to the Debye- Weiss 
theory, is increased. From Figs. 1-4 wc conclude that 
for hard core dipolar Yukawa fluids the MSA describes 
the linear and nonlinear dielectric permittivities with an 
adequate accuracy up to liquid number densities p*<0.8 



(note that freezing occurs at p* ~ 0.8; see Ref. [4t| for 
the related Stockmayer system), while the Langevin and 
the Debye- Weiss theory can be used only for p*< 0.1 and 
p*<0.2, respectively. 

Within the framework of our MSA theory and in agree- 
ment with our MC simulation data we have obtained a 
negative nonlinear dielectric permittivity for hard core 
dipolar Yukawa fluids. We note that using a cluster ex- 
pansion and the quadratic hypernetted chain (QHNC) 
approximation Martina and Stell 20[ obtained positive 
nonlinear dielectric permittivities for dipolar hard sphere 
fluids. They tried to explain the positive sign in terms 
of electrostriction (i.e., the dependence of p on E) which 
is included in their theory However, from an experimen- 
tal point of view this explanation is unlikely because the 
effect of electrostriction is so weak [i^] that it cannot 
overcompensate the contribution to Ae of molecular ori- 
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FIG. 4: Same as in Fig. 3 for T* = 2. Within the Debye- 
Weiss theory both eo and A* diverge at p* = 1.432. 



obtained: 

1) We have applied an extension of the mean spherical 
approximation (MSA) theory to determine the internal 
electric field dependence of the polarization for hard core 
dipolar Yukawa fluids. 

2) From the electric field dependence of the polarization 
analytical equations have been obtained for the linear 
and the nonlinear dielectric permittivity of these dipolar 
fluids. The predicted nonlinear dielectric permittivity is 
negative, which corresponds to the so-called normal sat- 
uration effect of dielectric media. 

3) Canonical Monte Carlo simulations have been carried 
out for the determination of the linear and nonlinear di- 
electric permittivity of dipolar Yukawa fluids. There is 
good agreement between the results from the MSA and 
the Monte Carlo simulation data for the reduced dipole 
moments (m*) 2 < 1 studied here (see Figs. 1-4). 

4) We have compared our theoretical results with the cor- 
responding Langevin and Debye-Weiss approximations. 
The observed trends are in agreement with the low den- 
sity behavior of the various approximations (see Eqs. 
© and I®). 

Our new theoretical approach provides a quantitatively 
reliable description for the nonlinear dielectric permittiv- 
ity. This raises the expectation that actual experimen- 
tal systems can be analyzed and understood along these 
lines. It will be also interesting to detect model systems 
which exhibit anomalous dielectric saturation, i.e., a pos- 
itive nonlinear dielectric permittivity. 



entations in the presence of an applied external field (see 
Ref. |3J|). Note that Eq. (33|) as well as Eq. |16]) arc also 
valid for dipolar hard spheres, which implies that within 
MSA the dispersion forces do not influence the dielectric 
properties. 
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